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Dynamic wetting on the flexible hydrophilic pillar- arrays is studied using large scale molecular dynamics 
simulations. For the first time, the combined effect of the surface topology, the intrinsic wettability and the 
elasticity of a solid on the wetting process is taken into consideration. The direction-dependent dynamics of 
both liquid and pillars, especially at the moving contact line (MCL), is revealed at atomic level. The flexible 
pillars accelerate the liquid when the liquid approaches, and pin the liquid when the liquid passes. The liquid 
deforms the pillars, resulting energy dissipation at the MCL. Scaling analysis is performed based on 
molecular kinetic theory and validated by our simulations. Our results may expand our knowledge of 
wetting on pillars and assisting the future design of active control of wetting in practical applications. 



Dynamic wetting on the hydrophilic pillar-arrays, where physics 1,2 , chemistry 3,4 , biology 5,6 , and materials 
science 7,8 intersect, is of great interest. Influenced by the combined effect of the surface topology and the 
intrinsic wettability of the solid, dynamic wetting on the pillars becomes rather complicated 9 . The topo- 
logy enhances wettability, making the substrate superhydrophilic 10 . As shown in Fig. l(a, c), a droplet on pillars is 
divided into two parts: the lower fringe penetrates into the spaces between pillars, whilst the upper bulk droplet 
spreads on the fringe. There are two important quantities to describe the pillar- arrays: surface roughness 
ro = 1 + 4dh/p 2 and density of the roughness (j) s = d 2 /p 2 , where d, h andp are the size, the height and the period 
of pillars, respectively. 

However, the pillars are rarely rigid, but usually flexible. In the last decade, the flexible pillar-arrays have been 
extensively used to study the properties of the cell-substrate and cell-cell interface 11 " 13 . These flexible pillars, especially 
with large slenderness ratio of hid or low Young's modulus 7, are sensitive to external forces. In Fig. l(b, d), 
the tip deflection of the flexible pillars near the moving contact line (MCL) caused by the interfacial tension is much 
larger than that away from the MCL. Therefore, to utilize and design the pillars in real applications, it is important to 
understand the dynamics and mechanisms of wetting on the flexible pillar- arrays. 

Dynamic wetting of a droplet on flexible pillar-arrays is far from being well understood. Ishino et al. studied the 
wicking process of a liquid film into pillars based on the Washburn law 14 . The MCL at microscale propagates in a 
zigzag way, which composes a smooth leading edge at macroscale 10 . For hydrophobic pillars, the surface rough- 
ness 15,16 , size of the droplet, line tension 17 , and the maximal three-phase contact line attainable along the actual 
droplet boundary 18 influence the pinning force and decide whether the MCL is sticky or slippery on the pillars 19 . 
For hydrophilic pillars, the topographic features 20,21 , the intrinsic wettability of the solid 22 " 24 , and the precursor 
film 25 affect the dynamics of wetting on the pillars and the shape of liquid spreading. Further, pre-deformed pillars 
would lead to uni-directional spreading 26 . However, as far as we are concerned, there is no study on the dynamic 
wetting of a droplet on the flexible pillars though it is widely used, especially in the field of biomedicine and lab- 
on-a-chip. 

In this paper, large scale molecular dynamics (MD) simulations are employed to explore the dynamics of a 
droplet on flexible hydrophilic pillars. Influenced by the combined effect of the surface topology, the intrinsic 
wettability and the elasticity of a solid, liquid superwets the pillars. The motions of both liquid and pillars, 
especially at the MCL, are revealed at atomic level. Flexible pillars accelerate the liquid when the liquid approaches 
and pin the liquid when the liquid passes. The spreading liquid deforms the pillars, resulting energy dissipation at 
the MCL. Based on the molecular kinetic theory (MKT), the scaling law for a droplet on flexible pillar-arrays is 
obtained and validated by our simulations. The spreading velocity is direction -dependent and is visualized by the 
flow pattern. Our results may help in understanding the dynamic wetting on flexible pillar- arrays and assisting 
the future design of active control of wetting in real applications. 
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(a) t= 0.0 ns (b)/=0.5ns (c)^l.Ons (d)/=1.5ns 



Figure 1 | Illustrations of a droplet wetting (a) on rigid pillar-arrays, and 
(b) on flexible pillar-arrays, (c-d) Experiments of a droplet wetting on 
flexible polydimethylsiloxane (PDMS) pillar-arrays (Y = 2.2 MPa). 

Results 

Large scale MD simulations implemented in LAMMPS 27 were car- 
ried out to explore the dynamic wetting of a droplet on a flexible 
hydrophilic pillars. The simulation domain is composed by flexible 
hydrophilic pillar- arrayed substrate (Y = 73.9 GPa) and a water 
droplet (radius R 0 = 3 nm). Topological parameters of the pillars 
are listed in Table 1. The droplet was modeled using the simple point 
charge/extended (SPC/E) water model with viscosity fi = 
0.729 mPa-s 28 , density p = 994 Kg/m 3 , surface energy y LV = 
0.0636 N/m 29 , which are close to those of real water. In order to 
accelerate the spreading, a bulk force F B = 1 nN was applied to the 
droplet by adding a constant force to every oxygen atoms along the 
negative z-direction. Considering both the governing effects and the 
computational efficiency, the value of F B is carefully tested and cho- 
sen not too large to overcome the effect of surface tension, as well as 
not too small to accelerate the simulations. The evaporated water 
molecules would quickly saturate the simulation box, so the evap- 
oration effect is ignored. 

When brought into contact with the hydrophilic pillars, a droplet 
collapses into two parts: the fringe penetrates into the spaces between 
pillars, and the bulk droplet spreads on the fringe. Meanwhile, the 
flexible pillars bend due to the spreading liquid. As shown in 
Fig. 2(b), the tip deflections labeled by the red arrows all point to 
the centre of the droplet, making an increase of the local roughness, 
and accelerating the liquid spreading. The interior corners formed by 
the pillars and the substrate provide additional driving forces to the 
liquid, making the velocity at the corner faster than the rest of the 
liquid 30 . From the top view of Fig. 2(c), we find that the initially 
circular droplet accommodate the arrangement of pillars and spreads 
in the shape of a square. Once the liquid reaches the pillars, the pillars 
are dragged towards the liquid. In Fig. 2(d), the fringe of the droplet 
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Figure 2 | Cross-section view and top view of a droplet on the flexible 
pillars (ro = 3.47). The cross section is labeled by the black dashed line. 
The red arrows are 8 times long of the tip deflections of the pillars. The blue 
balls are solid atoms. The liquid surface is shown in gray. 

advances in the space between the pillars and climbs up the pillars. 
The bulk droplet spreads on the base of the fringe, whilst the height of 
the bulk droplet decreases. Some of the pillars are dragged closer and 
interact with one another. The tip deflection near the MCL is 
approximately one order of magnitude larger than that away from 
the MCL. Especially when the liquid reaches flexible pillars, the 
deformation of these pillars is larger than those submerged in the 
liquid, which implies energy dissipation at the MCL 31 ' 32 . When the 
liquid advances to pass these flexible pillars, these pillars also deform 
to pin the liquid preventing its further spreading. 

We statistically obtained the evolution of the averaged radius of 
fringe .Ry with respect to time t, shown in Fig. 3(a-b). For the purpose 
of comparison, the spreading of a droplet on a smooth substrate is 
shown in black squares obeying a scaling law of R ~ t 117 (red dashed 
line), which is validated by theory, experiments and simulations 10 ' 33 . 
Fig. 3 exhibits two cases, the pillars with (a) h = 2.04 nm and (b) h = 
3.06 nm. In both cases, the spreading of a droplet on the pillars obeys 
approximately a scaling law of Rf ~ f 1/3 (green line) in the early stage. 
The spreading speed is found to increase with the decrease of <j> Si 
shown with the increase in darkness. Then, due to the energy dis- 
sipation in deforming flexible pillars, the velocity and the scaling 
exponent decrease. The number of pillars in contact with the liquid 
increases with Rj<j) s . Therefore, the pillar-arrays with larger <j> s dis- 
sipate more energy at the MCL, which causes an early decrease of the 
scaling exponent. 
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Figure 3 | Evolution of the fringe radius with respect to time, (a) h = 

2.04 nm and (b) h = 3.06 nm. The black squares represent the radius of a 
droplet spreading on a smooth substrate. The red dashed and green solid 
lines represent the scaling laws of R ~ t n and R ~ f 1/3 , respectively. 
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Discussion 

The dynamic wetting of a droplet on the flexible pillars is influenced 
by the combined effect of the surface topology, the intrinsic wett- 
ability and the elasticity of solid. To reveal the physical mechanisms, 
we adopt MKT to perform scaling analysis in the early stage of the 
dynamic process. 

The MKT model is a bottom-up model, which originates from the 
FrenkePVEyring 35 view of liquid transport as a stress -modified 
molecular rate process and is first proposed by Blake et al. to describe 
the MCL 36 . The motion of MCL is determined by the statistical 
dynamics of the molecules within the three-phase region. The energy 
dissipation is caused by the adsorption-desorption of the liquid 
molecules on the solid surface, and represented by the effective 
dynamic friction occurred at the MCL. In fact, when the liquid moves 
across a solid surface, the solid adsorbs and trends to immobilize the 
liquid molecules, while the liquid molecules desorb and tend to ad- 
vance. Considering the equilibrium liquid molecules in the vicinity of 
a solid, their advancing frequency k + and receding frequency rc~ are 
the same: k + = k~ = K 0 . k 0 is the equilibrium frequency and could be 
expressed in the term of the activation energy AG, K 0 = (k B T '/h) 
exp( — AG/NaIcbT), where k B , T, h and N A are the Boltzmann con- 
stant, absolute temperature, the Planck constant and the Avogadro 
number, respectively. The first term k B T/h is termed the fundamental 
frequency, which is a critical component and fundamental to any 
activated process in transition state theory. The second term 
exp( — AG/N A k B T) represents the influence of activation energy to 
the fluid molecules. In the case of liquid molecules adsorbing on a 
solid substrate, the activation energy AG is related to the work of 
adhesion between the solid and the liquid Wa and the spacing 
between surface sites A 32 ' 36 . The fluid molecules in the three-phase 
region are hindered not only by interactions with solid surface (adhe- 
sion), but also by viscous interactions with neighboring fluid mole- 
cules. So the activation energy AG = AG a + AG V , where AG fl = N A k 2 
Wa and AG V are the contributions from adhesion and viscosity, 
respectively. The relationship between the viscosity and AG V is 
fi=(h/v m ) exp(AG v /N A k B T) , where /n is the viscosity and v m the 
molecular flow volume, respectively. Hence, we could obtain 
K 0 = (k B T/fjLVm ) exp ( - k 2 Wa jk B T) . 

When a driving work per unit area w is applied to the liquid, 
the equilibrium is disturbed. The advancing frequency k + = 
k 0 exp(wk 2 /2k B T) and the receding frequency 7C~=K;oexp 

( — wk 2 /k B T) would no longer be in balance, and the contact line 
moves relative to the solid. The driving force tilts the potential surface 
and makes it more difficult for molecules to move in the receding 
direction. The resulting velocity U of the MCL could be expressed by 

,4. \ i k B Tk ( k 2 Wa\ wk 2 
U={k + -K~)k = 2 — — exp -— — sinh— — 1 
V ; fiVm V hT J 2k B T v ; 

When the driving work wk 2 is much less than the thermal energy k B T, 
the linear approximation sinh (w k 2 /2k B T) ~wk 2 /2k B T could be 
used. Hence, we obtain U = w/C= (wk 3 j fiv m ) exp( — k 2 Wa/k B T), 
where £ = k B T/K 0 k 3 is the coefficient of wetting-line friction. 

In the framework of MKT, the energy dissipation results from the 
friction at the MCL, measured by Wa, which has a link with the 
Young's modulus Y. A variety of solids can be modeled as LJ solid 37 , 
such as the face-centered-cubic (FCC) solid 38 used in our simula- 
tions, in which only the nonbonding interactions are considered. 
From the atomic aspect of view 39 , Wa is expressed by 
Wa~A S -L/dl~es-LG 6 s-LPsPL/dl> where A S _ L is the Hamaker 
constant between solid and liquid, p s and p L are the number density 
of the solid and the liquid, d L is the size of the liquid molecule, 
Sx-y = is the depth of the LJ potential well, 

Gx-y = (g X -x + g Y - y)/2 is the zero -crossing distance. Y of the LJ 
solid is Y~ A S -s/d 3 s ~ 8s-sG 6 s _ s p 2 s /d 3 s 39 , where d s is the lattice con- 
stant of the solid. Substituting parameters used in our simulations, 



we obtain Wa — 0.1 J/m 2 and Y ~ 100 GPa. Comparing their 
expressions, Wa monotonically increases with the increase of Y. 
However, there is no analytical solution of the relation between 
Wa and Y. Hence, we used the MD simulations to obtain their 
relation by changing s S -s and found a roughly linear increase of 
Wa with respect to Y, i.e. Wa ~ kY, shown in Fig. 4(a). 

The topology of the pillar- arrayed substrate plays a pivotal role in 
the wetting process. Since the evaporation of the liquid can be 
ignored in the MD simulations, the mass of the droplet remains 
constant, V 0 = V hu]k + y fringe , i.e. 4nR 3 0 /3 = nH b (3R 2 b + H 2 )/6 + 
n(l—(l> s )hR 2 , where R 0 , H b , R b , i^and <j) s are the initial radius of 
the droplet, the height of the bulk droplet, the radius of the bulk 
droplet, the radius of the fringe and the density of roughness (<j) s = 
d 2 /p 2 for pillars). We assume that otR b ~ Rf ~ R in the early stage, 
which is validated by experiments 10 ' 40 . Taking account of the lubrica- 
tion approximation 41 , i.e. H b < R b . We get 0 ~ H b /R b ~ 4R 3 0 /R 3 - 

3(\ — <j> s )h/R. Hence, a characteristic length L=^R? Q /h can be 
defined, which represents the maximum propagation radius of the 
fringe, to distinguish two extreme regimes 10 . When L^> R, the effect 
of surface roughness is neglected. The spreading of the droplet is only 
controlled by the interfacial energy and the bulk force, which is 
discussed in previous references 2 ' 33,42 and is not the focus in this 
paper. In our simulations, L <C R, i.e. relative small droplet, long 
or dense pillars, the effect of surface roughness dominates the 
dynamic wetting process. Not only the interfacial energy and the 
bulk force, but also the surface topology should be considered. 
Therefore 0~(\-$ s )h/R. 

In the dynamic wetting process, the driving work per unit area w 
consists of three parts: w = w y + w B — w Y . We discuss these three 
parts as follows: 

(1) w y is the driving work per unit area arising from interfacial 
tension. w y equals to the total change of interfacial energy when 
droplet spreads w y = (y sv — ysL)' ro ~yLv cos ft where y SL , y sy , 
y LV and 9 are the solid-liquid, solid-vapor and liquid-vapor 
interfacial energies and the instant contact angle, respectively. 
Taking account of the Young's equation (y sv — y SL = y LV cos0 0 , 
6 0 is the static contact angle) 43 , w y ~y LV (w cos 0 0 — cos 0). In 
the equilibrium state (w y = 0), for a rough surface (ro > 1) we 
obtain the Wenzel's relation: cos 9 = ro- cos 0 0 (0 < 0 0 ) 44 , for a 
smooth surface (ro = 1) we get back the Young's equation. 




0.0-1 ■ 1 ■ 1 ■ 1— 

0 100 200 300 
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Figure 4 | (a) The relationship between the work of adhesion Wa and the 
Young's modulus Yin the ordinary range of FCC LJ solids. The red line is a 
linear fit of the MD results, (b) Evolution of i^with respect to f for different 
Y. 
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Taking account of the lubrication approximation, w y ~y LV 6 2 
~y LV (\—(j) s ) 2 h 2 /R . In our simulations, w y ~ 10" 2 J/m 2 , and 
decreases with the expansion of the radius. 

(2) w B is the driving work per unit area arising from bulk force. 
Wb = FbH I Rj, where F B and H are the total bulk force and the 
height of the centre- of-mass, respectively. Set the surface of the 
substrate to be zero potential of the bulk force. Assuming the 
bulk droplet be a spherical cap, the height of the centre-of-mass 
of the bulk droplet equals h+ (2R 2 b H b + Hl)/(6R 2 b + 2H 2 h ) 
~ h + H b /3. We get w B ~ F B ( 1 - (j> s )h 2 /R* . In our simulations, 
w B ~ 10" 2 J/m 2 and mains a constant. 

(3) w Y is the resistant work per unit area arising from deformation 
of the pillars. w Y ~<j) s Yld 2 j h 3 p 2 , where I = d 4 /12 and d are the 
moment of inertia and the average deformation of the pillars, 
respectively. In our simulations, w Y ~ 10~ 3 J/m 2 in the early 
stage, one order of magnitude less than w y or w B . So we neglect 
Wy in the following derivations of scaling analysis in early stage, 
but the effect of deformation of the pillars will be considered in 
the last of this paper. 

From the above discussions, w = w y + w B — w Y ~ 10~ 2 J/m 2 , the 
driving work wl 2 is one order of magnitude less than the thermal 
energy k B T. Hence, the linear approximation could be adopted 
sinh(w2 2 /2A;BT) ~ wk 2 /lk^T. Substituting the above results into 
equation (1), we obtain 



dR 

U = — ~ exp 
dt r 



hT 



YLvii-lfh 2 



R 2 



(2) 



From equation (2), we get a dimensionless solution R/R 0 ~ 
(t/Tc) 1 ^ 3 , where the characteristic time T c = fiRl exp (A 3 Y / k^T^) / 
[y L y(l — (j) s ) 2 h 2 ] is controlled not only by the liquid properties 
(7lv> V an d Ro), but also by properties of the pillars (7, (j) s and h). 
When a droplet spreads on a smooth solid, the dimensionless solu- 
tion is R/Ro ~ (t/zj) 1 ^ , and the characteristic time is T c f = R 0 /v Ca 
= fiRo/y LV , which is only controlled by the liquid properties 10 ' 45 . 

Expanding the expression of R with respect to f, we obtain 



exp 



)?Y 



y LV (l-<fr s yk 



2^2 



(3) 



In the case of the same h and Y, (j) s determines the spreading 
velocity. The smaller (j) s is, the faster the liquid spreads on the pillars. 
Validated in Fig. 3 (the darkness increases with the decrease of </> 5 ), 
the spreading of the fringe increases with the decrease of <j) s . 

In the case of the same h and </> s , Y determines the spreading 
velocity. For a hydrophilic LJ solid, the smaller Y is, the faster the 
liquid spreads on the pillars. To validate this point, we simulated a 
droplet on the flexible pillars by only changing Y and keeping the 
other quantities constant. The topology in these MD simulations is 
the same as sample 2 in Table 1. And we obtained the evolution of Rf 
with respect to t for different Y in Fig. 4(b). Obviously, the spreading 
of the droplet increases with the decrease of Y. The increase of Y 
means an increase of Wa, making the pinning to the liquid molecules 
and the friction at the MCL be stronger. When Y is small, i.e. Y = 
36.93 or 73.90 GPa, the liquid spreading is smooth, obeying a scaling 
law of R ~ t m shown in Fig. 3(a). When Y is large, i.e. Y = 147.8, 
221.7 or 369.4 GPa, although the curves obey the scaling law of R ~ 
t m , they are not smooth. Large Y also implies large solid-liquid 
interaction. The pillars accelerate liquid when liquid approaches 
the pillars and also pin liquid when liquid passes the pillars, which 
makes fluctuation steps in the propagation curve. 

Fig. 5(a) shows a droplet propagating on a smooth hydrophilic 
surface. Fig. 5(b) shows the evolution of droplet radius in polar 
coordinate, in which the radial, the angular coordinate and the colour 
represent the time evolution, the direction and the spreading 
distance labeled by the below colour legend, respectively. In 
Fig. 5(b), point A (1.5 ns, 20°, orange) represents the liquid at 20° 




Figure 5 | (a) Snapshots of the wetting of a droplet on a smooth 
hydrophilic surface, (b) Evolution of droplet radius R in polar coordinate. 
The origin O is put on the centre of the droplet. The radial and the angular 
coordinate represent time evolution and the direction, respectively. The 
colour represents the spreading distance labeled by the below colour 
legend. The total time is 2.5 ns. (c) Evolution of R with respect to tat 45° 
labeled in (b). 

advances about 40-45 A in 1.5 ns. Line at 45° in Fig. 5(b) represents 
spreading radius of liquid at 45° with respect to time shown in 
Fig. 5(c), which obeys a scaling law of R ~ t in . The spreading of 
droplet on a smooth surface is nearly isotropic, i.e. independent of 
direction. 

The spreading on the pillars depends on the arrangement of the 
pillars, hence is direction-dependent. The flow patterns for samples 
1-4 are visualized correspondingly in Fig. 6(a-d). Because of the 
topographic features of the substrate, there exist fast and slow direc- 
tions of the liquid flow. However, the scaling law of R ~ t 1/3 is obeyed 
in all directions. Because the pillar arrays are symmetric, the flow 
pattern repeats every 45°. Comparing Fig. 6(a-d), we find the wett- 
ability gradually increase from Fig. 6(a) (more blue colour) to 
Fig. 6(d) (more red colour), which implies an increase of the overall 
velocity with the decrease of <j> s . In Fig. 6(a-d), there exist obvious 
boundaries distinguished by colors. The liquid advances faster in the 
region with more red colour and advances slower in the region with 
more blue colour. The boundaries of the fast and slow directions are 
sharp in Fig. 6(a), and gradually turn blurred in Fig. 6(d), which 
implies that the direction-dependence effect decreases with the 
decrease of roughness. 

In the last part, we discuss about the deformation of the pillars, 
which is neglected in the derivation of MKT in the early stage of 




Figure 6 | The direction-dependence of .R/with respect to f. Samples 1-4 
(a-d) with flexible pillars, and (e-h) with rigid pillars. The origin O is put 
on the centre of the droplet. The radial and the angular coordinate 
represent time evolution and the direction, respectively. The colour 
represents the spreading distance labeled by the right colour legend. The 
total time is 2.5 ns. 
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Figure 7 | The averaged displacement of and force on each layer of the pillars for (a) ro = 2.65, (b) ro = 1.89, (c) ro = 2.78 and (d) ro = 2.04. The gray 
and orange circles represent the initial and final positions, respectively. 



spreading. However, the spreading velocity of the liquid is influenced 
by the deformation of the pillars 46 . To compare with the propagation 
of samples 1-4 with flexible pillars, we set the pillars to be rigid in the 
MD simulations and obtain Fig. 6(e-h). The overall velocity of 
Fig. 6(a-d) (more blue colour) is less than that of Fig. 6(e-h) (more 
red colour), because deformations of these pillars dissipate energy. 
With the increase of </> 5 , the energy dissipation by the flexible pillars 
increases. Although the stiffness is changed, the fast and slow flow 
directions are nearly the same. The direction-dependence of the 
velocity is not changed. 

The forces applied on and the corresponding deformation of each 
layer of the pillars varied with the roughness, shown in Fig. 7. To 
eliminate the thermal fluctuations, the forces and displacement of 
atoms were averaged over all the pillars and the time. In the case of 
large roughness (Fig. 7(a, c)), pillars are close to each other. The fringe 
reaches the upper and lower parts of the pillars almost at the same 
time, resulting in an approximately uniform distribution of the forces. 
The pillars bend into a quadratic curve. In the case of small roughness 
(Fig. 7(b, d)), pillars are far from each other. The fringe firstly reaches 
the lower part, then the upper part after a short time, resulting in the 
forces exerted on the lower layers be a little larger than those on the 
upper layers. The pillars bend into a sigmoidal curve. 

In summary, the dynamics of a droplet on the flexible hydrophilic 
pillars is explored using large scale MD simulations. For the first 
time, the combined effect of the surface topology, the intrinsic wett- 
ability and the elasticity of a solid on the dynamic wetting is taken 
into consideration. The atomic details of both the liquid and the 
pillars during the superwetting, especially at the MCL, are revealed. 
The flexible pillars accelerate the liquid when the liquid approaches, 
while also pin the liquid when the liquid passes. The direction- 
dependent spreading of the liquid deforms the pillars, resulting in 
energy dissipation at the MCL. Our results may help to understand 
the dynamic wetting of the droplet on flexible pillar- arrayed sub- 
strate and assist the future design of active control of wetting in 
practical applications. 

Methods 

Molecular dynamics simulations. Large-scale MD simulations implemented in 
LAMMPS 27 were carried out to explore the dynamic wetting process of a droplet on 
flexible hydrophilic pillar-arrays at the atomic level. The total potential energy Eg 
between two atoms i and j separated by is the sum of Lennard- Jones (LJ) potential 
energy and the Coulombic pairwise interaction 



= 4e« 



(4) 



where e is the depth of the potential well, a is the zero-crossing distance for the 
potential, k e = 8.988 X 109 N-m 2 /C 2 is the Coulomb constant and q is the charge on 



atom. The cut-off length for LJ potential is 1.5 nm, while that for Coulombic 
interaction is 1.0 nm. 

The NVT ensemble (constant number of atoms, volume and temperature) was 
used. The Nose-Hoover thermostat, acting as a heat reservoir, with a time-step of 1 fs 
was employed to regulate the temperature at 350 K. The simulation domain is about 
15 X 15 X 20 nm 3 according to the arrangement of pillars. 

The water molecules were simulated using the simple point charge/extended (SPC/ 
E) water model: the oxygen atoms were modeled as charged LJ particles (o"o-o = 
0.3166 nm, £ 0 -o = 0.650 kj/mol, q Q = —0.8476 e), while the hydrogen atoms were 
charged but without considering LJ interactions between them (cr H _ H = 0.0 nm, 
s H -h = 0.0 kj/mol, q H = 0.4238 e). The solid atoms were modeled as uncharged LJ 
particles with (T S - S = 0.2637 nm and s S -s = 42.5723 kj/mol, which were chosen to 
be hydrophilic to the SPC/E water. The values of o and s between them were calcu- 
lated according to the Lorentz-Berthelot rule: o X -y = (cx-x + Cy-y)/2 and 
b X -y = y/Gx-xGY-Y> where the subscripts X and Y represent the types of atoms. So 
(t 0 -h= 0.1583 nm,<7 s _ H = 0.1319 nm,<7 s _ 0 = 0.2902 nm and s 0 -h = 0.0 kj/mol, 
e S -h = 0.0 kj/mol, e S -o = 5.2604 kj/mol. The LJ interactions between hydrogen and 
oxygen atoms, as well as between hydrogen and solid atoms were ignored. The cut-off 
length for LJ potential is 1.5 nm, while that for Coulombic interaction is 1.0 nm. 

The SPC/E model was established by adding an average polarization correction to 
the SPC model, with a modified value of q Q (charge of oxygen atom) and q H (charge of 
hydrogen atom). Compared with SPC, transferable intermolecular potential 3 point 
(TIP3P), TIP4P, TIP5P water model, SPC/E model predicted a better density (p = 
994 Kg/m 3 ), viscosity (/i = 0.729 mPa- s) 28 and surface tension (y LV = 0.0636 N/m) 29 , 
which are close to those of the real bulk water. 
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